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A general method to derive the master equations for extremal models is established. These 
systems are shown to develop a peculiar kind of correlations between elements related to the char- 
0^ ' acterization of extremal dynamics as an information process. 
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I. INTRODUCTION 



In this article we would like to establish the general method and the proper way to derive the master equations 
governing the Bak-Sneppen type models of extremal dynamics El^. Apart from the relevance as rough models of 
biological evolution, the interest of these systems is that they show at equilibrium a self organized critical state, with 
■ power law correlations in both time and space (if any 'space' relation between elements is defined, as in the case 
of the nearest neighbor model) and the phenomena of avalanches and punctuated equilibrium. The importance of 
the extremal dynamics in determining the self organized critical properties of several models (another example is the 
Invasion Percolation) has been discussed in H], 

The systems under consideration are composed of a set of, say, L elements characterized by the value of a variable, 
usually a one dimensional variable taken in the unit interval (in the evolutionary interpretation it parametrizes the 
facility to mutate). Let's denote it xa where A = l,...,L is the element label and xa € [0, 1]. The characteristic 
feature of the dynamics of these systems is the selection, at each update, of the element which has the maximum 
value (the "maximum element") and the reassignment to it of a new value randomly chosen in the unit interval. This 
kind of rule is referred to as extremal dynamics. In addition new random values are assigned to a certain number of 
other elements which, in the more interesting case, have a neighborhood relation with the maximum element (nearest 
neighbor models, n.n. models in the following). Another possibility is to choose them at random (random models), in 
which case the model is analytically solvable via a master equation of the diffusion type for the probability P\(n) to 
have n elements with values less than a given threshold A |Q. Whereas this latter model has the advantage of giving 
^ . exact results, it is not very rich, lacking any spatial structure and in some sense not being an interacting system, 
as the maximum element affects the other elements just by triggering the change of their values. To write a master 
equation for the n.n. model it is not anymore possible to work with the statistics of the variables x disregarding 
the precise assignments to the elements, as it is the case in the random model. Instead we have to cope with the 
full probability density p(x), x — (xi, . . . ,xl), that is with probabilities defined in the hypercube [0, 1] L . To our 
00 knowledge the master equation for the n.n. model has not yet been derived. 

0^ ■ Now an interesting point is the characterization of extremal dynamics as an information process: at each time 
step the choice of the maximum element takes into account the states of all the elements (or at least a large number 
of them: in Invasion Percolation the perimeter elements) giving information about them. As we'll see below, this 
global information creates correlations even between elements that remain inactive. We'll refer to them as information 
correlations. As pointed out in Q the role of fluctuations is essential to determine the properties of the self organized 
critical state of the Bak-Sneppen models, in particular the avalanche phenomenon and the distribution of the maximum 
element: in fact the relevant quantity, the number of active elements above self organized threshold, is of the same 
' order of magnitude as its fluctuations. This limits severely a mean field theory treatment where correlations between 
^ . elements are neglected, and that even in the case of the random model where they arise only as information correlations. 

X: 

II. HOW CORRELATIONS BETWEEN ELEMENTS ARE GENERATED 

The key feature of our method is to decompose the hypercube in a partition where the regions correspond to a 
given ordering of the element values. Let's start by introducing some definitions and notations. Iai,...,A„ will denote 
the domain of the hypercube where the values of the n elements A\ , . . . , A n are decreasingly ordered: 

I Al ,...,A n = {xe[0 1 l} L / x Al >..->x An >x B VB=£A 1 ,...,A n } (1) 
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and A.a 1 A„(x) will denote the characteristic function which has this domain as support: 

a ^.-.A»W-| 0) othcrwise . 

First we have to determine how the choice of the maximum element as the active one affects the probability distribution. 
If at a given time step the probability distribution is p{x) , the knowledge that the maximum element is the element 
T changes p in the following way: 

p(x) -» P (x\I T ) = NA T (x)p(x). 

The normalization factor N, which depends in general on T, is the inverse of the probability to have T as the maximum 
element and N — L if the function p is symmetric. If we now replace its value with a new one randomly chosen in 
the unit interval and with probability distribution ip(x), the updated distribution function becomes: 

p(x) p(x;T) — <p{xt) I dxT p(x\It) = N(p(xr) / dx T p(x). (2) 

JO Jmax{x A ;A^T} 

In the following we'll consider the case <p(x) = 1. Now we want to show that this new distribution represents correlated 
elements even if p(x) does not, that is the knowledge of which is the maximum element creates correlations between 
the other elements. To see this let's consider L = 3 and p(x,y, z) = 1. If we are now given the information that the 
first element, corresponding to x, is the biggest one, then the probability distribution restricted to y and z is: 

p(y, z)= f dx p(x, y, z\h) = N [6{y - z)(l - y) + 8(z - y)(l - z)} 
Jo 



and 

p{z) = I dy p{y,z) = N 
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Then the conditional probability p(y\z) = p(y,z)/p(z) depends on z, that is y and z are no more independent: an 
information on one of the two, let's say z, gives information on the other. To understand this suppose z = £; this 
means x > £ and, because y < x, we have p(y\z = £) = const, for y < £. On the other hand p(y\z = £) has to tend to 
for y — ► 1: if y = 1 — e then x > 1 — e and the probability of this last condition tends to zero for e — * 0. This result 
contrasts with what is said in , where it is erroneously assumed that the element variables are independent after 

the application of the extremal rule, and that the probability distribution maintains a factorized form. As we already 
pointed out the origin of these correlations comes from an information process and not from an action connecting the 
elements. This last mechanism is only present in the n.n. model where it is responsible of spatial correlations. It 
has been shown in B that in the random model mean field and exact results (supplied by simulations) have strong 
discrepancies and we see now that they can be attributed entirely to these information correlations. 



III. THE SIMPLEST MODEL 



To show the technique used to derive and cast in a convenient form the master equations of the Bak-Sneppen 
models, we examine in this section the quite simple case where the maximum element is the only one to be updated. 
The equation we are looking for is that which gives the probability distribution at time step i + 1, pi+i(x), in terms 
of that at time i, Pi(x). If Pi(T) is the probability to have element T as the maximum at time i and noting that it is 
the inverse of the normalization factor N in equation (0), the master equation for the simplest model is: 

l „i 

Pi+i{x) = ^2 m*( t )k(^; t ) = X! / dXT Pi ^- 

T=l j* J max{xA;A^T} 

Our strategy in investigating this equation is now to decompose the hypercube [0, 1] L in regions of the type ([!]) such 
that the probability distribution maintains an analytic form in each cell. Of course the finest partition where we take 
as cells the regions Ia 1 ....,a l where L is the total number of elements does the job for every model. In the very simple 
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case under consideration it turns however out that the partition in cells Ia,b is sufficient. Then we decompose the 
probability distribution in the following way: 



p(x) = J2A AtB (x)p {A ' B \x) 



A,B 

where the functions p( A ' B ^ can be extracted from p by using Aa,b as a projector 

p^ = Aa,b P 

and we want to express the functions p[ A ' B ^ in terms of the p^ C,D ^'s. It turns out that they remain analytic after the 
probability transition i — ► i + 1. To manipulate these " pieces" in the master equation where integrals are taken over 
domains of integration that cross the boundaries of the partition, we need to take into consideration other kinds of 
regions: 



^ x n = e [0, l] L I x Al > ■ ■ ■ > x An >x B \/B^D,A u ..., A n } 



and, as before: 



^M,...,A n {x) = 



1, iix&l 
0, otherwise. 



D 

Ai,...,A r 



By means of the rules: 



AtA a ,b = St,aAa,l 



Ad,eA]^ — Ad,e(<>b,d + Sb,e5t,d) 
we have, for a symmetric function p(x) and N = L: 

P (x\I T ) = LA T (x)p(x) =lJ2&tA a ,bP {A ' B) = lJ^ &t,bP (T ' B) 



A,B 



then 



and finally 



p(x; T)= f dx T p(x\I T ) =LJ2&1 [ dx T P (T ' B) 



1 f 1 

p i+ i(x) = -^2pi(x;T) =^A| / dx T p 

T T B Xb 



(T,B) 



( D > E ) A A 



£ f dxr p7' D) + f dx D pf^ 

'/J n> X D •> XE 



T(^D) JXD 



which is the desired result, that is an analytic form for p( D ' E "> in its domain of definition Id,e- If po = 1 this recursion 
equation is easily solved to give: 

p?> e \x) = ( ^ L 1 + L)! [{l - i)(i - XD y + (i - XE y] . 

The function Pi(x) encodes clearly the most detailed information on the system at time i. We can calculate from it 
for example the one-element probability distribution. Choosing the element 1: 



Pi(xi) 



/ dx 2 / 

Jo Jo 



dx 3 --- I dx L Pi(x) 
o 
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rxi rx^ rX2 

= {L - 1) / dx 2 / dx 3 ■ ■ ■ / rfx L p< 1,2) 
Jo Jo Jo 

+(£ - 1) / r dx 3 ■ ■ ■ f 1 dx L P f- i] 

JX! JO JO 

pi rx 2 rx 3 rx 3 

+{L-l)(L-2) dx 2 dx 3 dx 4 --- dx L p\ 

Jo Jx t Jo Jo 



The first term is the probability distribution (to be normalized) of the maximum element, the second is that of the 
next smaller one, and the third term the probability distribution of any of the others (in this case they are all equally 
distributed). By substituting the expressions we found for the p\ D,E ^ the result for the mean value is: 

/ /.nv f ,l 1 1 L-l 

= / d x x lPl = — + - ■ — — 
J[o,i] i lL 1 L + i 

which, apart from the term 1/(2L) coming from the contribution of the value replaced at each update, goes to zero 
as i~ x . 



IV. MASTER EQUATION FOR THE N.N. MODEL 

We move to the derivation of the master equation in the case of a simplified version of the n.n. model, where 
besides the maximum element only one of the two adjacent ones is updated, chosen with equal probability to preserve 
left-right symmetry. The partitioning into regions Ia.b is not anymore sufficient, in the sense that at each update the 
regions of analyticity of the probability distribution break up into finer and finer regions. Then we have to decompose 
into domains 1^ where for brevity we use the vector notation A = {A\, . . ., Al). To this vector we associate the 
permutation <JA{k) = Ak- Let's introduce some notations, operations on vectors A and rules. If we take away the 
element Ak from the vector A the resulting L — l dimensional vector is denoted as 

A(A k ) = {A u ...,A k - U A k+1 ,...,A L ) 

and if we take away element T and move it in the k position: 

A(T:k) = (A(T)i, . . . , A(T) k _i,T, A(T) k , A(T) k+ i, . . .). 

Some useful identities are: 

A T A A = S T,Ai A A(T:l) 



where F(B) is any expression involving vector B. Let's start by rewriting the result of the removal of the maximum 
element in terms of the finest decomposition that is p B (x; T) in terms of p B (x). Given that: 

p{x\I T ) = LA T (x)p(x) = L^ &a(T:1)P MT:1) 

A(T) 



we have 



p(x; T) = f dx T p(x\I T ) =LJ2 ^ T a(t) I dx ? P X(T:1 

J J X r™. 



A(T) • 4 < T >i 
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and then 

pi pi 



p S (x;T) =LJ2 %(T),A» / d *T P A{T - A) =L[ dx T p 8 ^. 



A(T) 



-A(T) 1 



The next step is to find the change in the p(x) after having removed an arbitrary element D. The new probability 
distribution is denoted p(x; [D]): 

p(x; [D]) =J2[ dx ° A S P 8 W = E A S(D) dx ° 

B B J ''BW 

where the operator 7rg on an element is given by ir B {Bk) — Bk+ n and by definition x^-i, B \ = 1, x w i ( Bl ^ = \/B. 
Projecting 



P A 



f X _1 (- D ) ~* 

{x;[D\) = ^ S A{D),S{D) j B ° dx D p B (x) 



E S A(D),B(D) 



l{D)k - 1 dx D p^Hx). 



^ , A ^ k y dx D P A( - D:k \x) 



E 

k " ^(D) 



Denoting the process of updating successively the maximum element and another one by p A (x;T, [D]) we have for 
the corresponding master equation: 



P?+iW = E T , D ^(T, [D])p A (x; T, [£>]) 



a 

k ^Ain) 



where in the case of the n.n. model the hat indicates sum over n.n., fii(T, [D]) = (2i) _1 and a = 1/2. This is a 
mapping for the set of L\ functions p A (xi, . . . , Xl). In the case of a system on a ring, due to rotational and reflection 
symmetries, there are only (L — l)!/2 independent functions. In fact if A and B are related by a translation (a cyclic 
permutation) and/or a reflection we have: p A (x) = p B (y) with XA k = VB k - The number of functions in the game can 
be further reduced, as we'll see in the next section. 

V. THE CASE L=4 

To test the correctness of the last master equation we'll apply it to the case L — 4 which is the minimum size able to 
generate information correlations. We choose as independent the functions pC 1 ' 2 ' 3 ' 4 ), pf 1 * 2 . 4 * 3 ) and pC 1 ' 3 ' 2 ' 4 ), Moreover, 
these three functions can be expressed in terms of a single two variables function, g(x,y), in the following way: 

p( 1 ' 2 > 3A \xi,x 2 ,x 3 ,x 4 ) = [g(xi,x 2 ) + g(x 2 ,x 3 ) + g{x 3l Xi) +g(xx,Xi)]/4 

p {1 ' 2 ' 4 -' 3 \xi,x 2 ,x 3 ,x i ) = \g(xi,x 2 ) + g{x 2 ,x 3 ) +g{x Al x 3 ) + g{x 1 ,x i )]/A 

p (1 ' 3 ' 2A ' ) (xi,x 2 ,x 3 ,x i ) = [g{xi,x 2 ) + g(x 3 ,x 2 ) + g{x 3 ,x A ) + g(xi, ar 4 )]/4. 
In terms of g the master equation (|^) becomes: 

g'{x,y) = - [(1 -x 2 )g(x,y) +g x (x,x) +xg x (x,y) + (1 - x)g x (y,y)+ 
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(1 - x)g y {x) + (1 - x)g y (y) + 2g xy (x) + 2G(x)} 

where we have defined 

9x(x,y) = / d£ g(£,y) 



g y (x) = dr) g{x,rj) 
Jo 

9xy{x) = dr) g x {x,ri) 



G(x) = / d£ g x (£,S). 



The two variables function g can be directly interpreted as a probability distribution: g(xA,XB) is the probability 
distribution of the values of elements A and B with the condition xa > xb after the update of the other two. This 
allows us to express the one-element probability distribution (without the trivial contribution coming from the two 
updated elements) as: 

f[x)= I dr)g{x,rj)+l d£ g(£, x) = g y {x) + g x {x, x). 

JO Jx 

Starting with the uniform initial condition, go(a,b) — 1, the master equation can be iterated exactly. The resulting 
mean value m, = f Q dx xfi(x) is in perfect agreement with simulations. 



VI. CONCLUSIONS 



Our main point was to show the presence and the origin of correlations related to the peculiar mechanism of extremal 
dynamics and to develop a formalism to treat them exactly in the case of Bak-Sneppen models. This leads in general 
to a linear mapping of L\ functions in L variables, the probability distributions of the element values subjected to one 
of the L\ possible orderings. This mapping can be further simplified by taking into consideration the symmetries of 
the dynamics as it was shown in the case of a four elements system where the 24 functions in 4 variables are reduced 
to only one in two variables. The hope is to reduce the generic finite size system to something tractable, at least as a 
starting point for approximations. 
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